clc;
clear;

idx = 1;
switch idx
    case 1
        V = [0 0 0;
            1 0 0;
            0 1 0;
            0 0 1];
        T = [1 2 3 4];
    case 2
        V = [0 0 0;
            1 0 0;
            0 1 0;
            1 1 0;
            0 0 1;
            1 0 1;
            0 1 1;
            1 1 1];
        T = [1 2 3 5;
            2 3 5 8;
            3 2 4 8;
            2 5 6 8;
            5 3 7 8];
    case 3
        D = [0, 1, 0, 1, 0, 1];
        N = 3;
        [V, T] = process_cube(D, N, N, N);
    case 4
        V = [0 0 0;
            1 0 0;
            0 1 0;
            1 1 0;
            0 0 1;
            1 0 1;
            0 1 1;
            1 1 1];
        T = [1 2 4 8;
            1 3 4 8;
            1 2 6 8;
            1 3 7 8;
            1 5 6 8;
            1 5 7 8];
    case 5
        V = [0 0 0;
            1 0 0;
            0 1 0;
            1 1 0;
            0 0 1;
            1 0 1;
            0 1 1;
            1 1 1;
            0.5, 0.5, 0.5];
        T = [1 2 4 9;
            1 2 6 9;
            1 5 6 9;
            1 3 4 9;
            1 3 7 9;
            1 5 7 9;
            2 4 6 9;
            3 4 7 9;
            4 6 8 9;
            4 7 8 9;
            5 6 7 9;
            6 7 8 9];
end

% showTetra(V, T)
pde.domain = [0, 1, 0, 1, 0, 1];
pde.bdary = [0, 0, 0, 0, 0, 0];
% mesh = process_mesh(pde, 1/2);
% fs = process_fespace(mesh, "P2");
% process_showmesh(mesh, fs, [1, 1, 0]);


% N = 2;
% [V, T] = process_cube(pde.domain, N, N, N);
% showTetra(V, T);
% [Tri, bdf, elem_tri, edge, elem_ed] = process_face_edge(V, T);

t = [4 3;
    3 2;
    2 4;
    1 3;
    4 1;
    2 1];
edge = reshape(T(:,t), 6, 2);
midpoint = 0.5*(V(edge(:,1),:) + V(edge(:,2),:));
Vb = [V; midpoint];
Tb = 1:10;
showTetra(Vb, Tb);
Tb = [1 5 6 7;
    1 5 6 8;
    1 5 7 9;
    1 6 7 10;
    2 6 10 7;
    3 8 6 5;
    4 5 9 7];
showTetra(Vb, Tb);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% process_face
function [Tri, bdf, Tetra_tri, Edge, Tetra_ed] = process_face_edge(V, T)
Ne = size(T,1);
tri_all = zeros(Ne*4,4);
t1 = [4 3 2;
    1 3 4;
    4 2 1;
    1 2 3];
edge_all = zeros(6*Ne,3);
t2 = [4 3;
    3 2;
    2 4;
    1 3;
    4 1;
    2 1];
k1 = 1;
k2 = 1;
for e = 1:Ne
    tri_all(k1:k1+3,:) = [e*ones(4,1), reshape(T(e,t1), 4, 3)];
    edge_all(k2:k2+5,:) = [e*ones(6,1), reshape(T(e,t2), 6, 2)];
    k1 = k1+4;
    k2 = k2+6;
end
[~, ib, ic] = unique(sort(tri_all(:,2:4),2), "rows", "stable");
k = accumarray(ic,1);
Tri = tri_all(ib,:);
bdf = find(k==1);
Tetra_tri = reshape(ic, 4, [])';
Hexa_tri = reshape(ic, 24, [])';
Nhex = Ne/6;
tri_face = zeros(6,Nhex+1);
t3 = [7 10;
    8 11;
    14 17;
    4 18;
    13 16;
    2 5];
tri_face(:,1) = (1:6)';
k = 2;
for e = 1:Nhex
    temp = unique(Hexa_tri(e,:), "stable");
    tri_face(:,k:k+1) = reshape(temp(t3), 6, 2);
    k = k+2;
end
k = ismember(tri_face(:,2:end), bdf);
temp = zeros(size(tri_face));
temp(k==1) = tri_face(k==1);
for e = 1:6
    fi = temp(e,temp(e,:)>0);
    Tri(fi,5) = e;
end
bdf = Tri(Tri(:,5)>0,:);
[~, ib, ic] = unique(sort(edge_all(:,2:3),2), "rows", "stable");
k = accumarray(ic,1);
Edge = edge_all(ib,:);
bde = find(k<=3);
Tetra_ed = reshape(ic, 6, [])';
end

%% process_refine
function [new_node, new_elem] = process_refine(node, elem, edge, elem_ed, basis_type)
Np = size(node,1);
switch basis_type
    case "P1"
        new_node = node;
        new_elem = elem;
    case "P2"
        midpoint = 0.5*(node(edge(:,1),:) + node(edge(:,2),:));
        new_node = [node; midpoint];
        e1 = abs(elem_ed)+Np;
        new_elem = [elem, e1];
    otherwise
        error("Invalid basis type.");
end
end

%% process_cube
function [Vectex, Tetra] = process_cube(domain, Nx, Ny, Nz, varargin)
[X, Y, Z] = ndgrid(linspace(domain(1), domain(2), Nx+1), linspace(domain(3), domain(4), Ny+1), linspace(domain(5), domain(6), Nz+1));
Vectex = [X(:), Y(:), Z(:)];
Npxoy = (Nx+1)*(Ny+1);
idxp = [0 1 Nx+1 Nx+2 Npxoy Npxoy+1 Npxoy+Nx+1 Npxoy+Nx+2];
[Hx, Hy, Hz] = ndgrid(1:Nx, 0:Ny-1, 0:Nz-1);
idx = reshape((Hx + Hy*(Nx+1) + Hz*Npxoy), [], 1);
Hexa = repmat(idxp, size(idx,1), 1)+idx;
if isempty(varargin)
    Tetra = zeros(6*Nx*Ny*Nz, 4);
    t = [5 1 7 8;
        1 5 6 8;
        2 1 6 8;
        1 2 4 8;
        3 1 4 8;
        1 3 7 8];
    k = 1;
    for i = 1:size(Hexa,1)
        Tetra(k:k+5,:) = reshape(Hexa(i,t), 6, 4);
        k = k+6;
    end
elseif varargin{1} == 5
    Tetra = zeros(5*Nx*Ny*Nz, 4);
    k = 1;
    layer = 1;
    flag_flip = 1;
    for i = 1:size(Hexa,1)
        row = mod(mod((i-1),Nx),Nx*Ny)+1;
        col = mod(floor((i-1)/Nx),Ny)+1;
        if mod(row+col,2) == 0
            flag = flag_flip;
        else
            flag = abs(flag_flip-1);
        end
        if flag == 1
            t = [1 2 3 5;
                2 3 5 8;
                3 2 4 8;
                2 5 6 8;
                5 3 7 8];
        else
            t = [1 5 6 7;
                4 1 6 7;
                1 2 4 6;
                3 1 4 7;
                6 4 7 8];
        end
        Tetra(k:k+4,:) = reshape(Hexa(i,t), 5, 4);
        k = k+5;
        if mod(i,Nx*Ny) == 0
            flag_flip = abs(flag_flip-1);
            layer = layer+1;
        end
    end
end
end

%% showTetra
function showTetra(V, T)
Np = size(V,1);
Ne = size(T,1);
figure("WindowStyle", "docked");
tetramesh(T, V, "FaceColor", "none");
hold on;
% node
plot3(V(:,1), V(:,2), V(:,3), "k.", "MarkerSize", 10);
text(V(:,1), V(:,2), V(:,3), num2str((1:Np)'), "FontSize", 10, "FontWeight", "bold", "Color", "k", "HorizontalAlignment", "right", "VerticalAlignment", "bottom");
% elem
centroid = zeros(Ne,3);
for e = 1:Ne
    E = V(T(e,:),:);
    centroid(e,:) = mean(E,1);
end
text(centroid(:,1), centroid(:,2), centroid(:,3), num2str((1:Ne)'), "FontSize", 10, "FontWeight", "bold", "Color", "r", "HorizontalAlignment", "right", "VerticalAlignment", "bottom");
hold off;
xlabel("X");
ylabel("Y");
zlabel("Z");
view(120, 30);
end
